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Abstract 

Efficient fourth order symplectic integrators are proposed for numerical integration of separable Hamiltonian systems H(p,q) = T{p) + 
V(q). Symmetric splitting coefficients with five to nine stages are obtained by higher order decomposition of the simple harmonic oscillator. 
The performance of the methods is evaluated for various Hamiltonian systems: Integration errors are compared to those of acclaimed integrators 
composed by S. Blanes et al. (2013), W. Kalian et al. (1999) and H. Yoshida (1990). Numerical tests indicate that the integrators obtained in this 
paper perform significantly better than previous integrators for common Hamiltonian systems. 
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I. Introduction II. Fundamental theorems 


Explicit symplectic integrators (Sis, singular SI) are numer¬ 
ical integration schemes for separable Hamiltonian systems. 
Sis are widely used in fields of celestial mechanics, accelerator 
physics, molecular dynamics, and quantum chemistry due to 
their structure-preserving properties and simple implementa¬ 
tion. 

Numerous Sis with increased accuracy and stability have 
been published since the first discovery of the fourth order 
SI [1]. The perhaps most acclaimed splitting composition is 
derived by H. Yoshida who proved that Sis of arbitrarily high 
even order exist and can be constructed using Lie algebra 
[2], Other more recent studies on the topic of high stability 
Sis include force gradient schemes [3] and symplectic correc¬ 
tors [4] [5], The integration order is often associated with the 
efficiency of numerical integration schemes. Nth order inte¬ 
grators have their integration error reduced by a factor k N 
when the time step is reduced by a factor k, making high or¬ 
der integrators highly favorable for small time steps. However, 
higher order error terms cannot be neglected for large time 
steps. Symplectic correctors have been proposed for eliminat¬ 
ing high order errors of Sis, although these methods are only 
competitive for specific systems. The force-gradient composi¬ 
tion has been proposed as a more efficient splitting scheme for 
general separable Hamiltonian systems but requires evalua¬ 
tion of the force gradient [3]. The force-gradient compositions 
benefit from all-positive splitting coefficients meaning that 
only forward steps are taken in phase space and are therefore 
referred to as a forward Sis. Non-gradient splitting methods 
of orders higher than two always involve negative coefficients 
[5]. To construct more efficient high order Sis it has been pro¬ 
posed that the absolute value of negative coefficients and the 
sum of positive coefficients should be minimized [5]. 

This paper presents a series of efficient fourth order Sis 
which are ideal for numerical integration of particle systems 
with quadratic kinetic energy. Additionally, some of the pre¬ 
sented integrators exhibit higher order accuracy for systems 
that resemble the harmonic oscillator. The presented splitting 
coefficients generally have small absolute values and are there¬ 
fore referred to as near-forward Sis. 
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Consider the general separable Hamiltonian H(q,p) = 
V{q) + T(p) with p = —dH/dq and q = dH/dp where q and p 
are canonical coordinates conventionally labeled position and 
momentum respectively. By defining z = ( q,p ), the time evo¬ 
lution of z can be expressed by [2]: 

z( T ) = e T M+ B ) z (o) ( 1 ) 

where A and B are non-commutative operators associated 
with T and V, and r is the time step. A set of real numbers 
exist (cj,C 2 , ...,C]t) and (di,d2, — ,*4), referred to as splitting co¬ 
efficients, that satisfy X^ =1 c,- = di — 1 such that the time 
evolution of z can be approximated by a product of exponen¬ 
tial functions [2]: 

e r(A+B) = YY^rA^rB + Q ^N+U (2) 

where o( t n+1 ) denotes an error on the N + 1 order term. It 
is common practice to expand the exponential terms by power 
series with respect to r. Lie algebra can be applied to de¬ 
termine the splitting coefficients which satisfy Eq. 2 exactly 
up to a given order. The time evolution of z (0) —> z (t) or 
(q 0 ,p 0 ) -A- (q k/ p k ) is explicitly computable as [2]: 

Vi = Pi-i ~ zd i% (<?;) 

for i = 1,2, ...,k. Note that the intermediate steps ( qj,Pj)f = 
1,2, ...,k — 1, referred to as substeps, do not represent the time 
evolution of (qo,po). Sis are time reversible if the symplectic 
coefficients are symmetric [2]. Time reversibility implies that 
z (0) —> z (t) can be exactly reverted by z (0) z (— r). This 
property is embraced throughout this paper. The symmetric 
composition in the form of Eq. 3 with d k = 0 is referred to as 
ABA [5]. If the roles of A and B in Eq. 2 are switched the time 
evolution is explicitly computable as: 

ar / \ W 

<\i = Vi-l + TC, ^ ( : pi ) 

for i = 1,2,..., k. The symmetric composition in the form of Eq. 
4 with c k = 0 is referred to as BAB. Table 1 summarizes the 
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symmetric structures of the ABA- and BAB-composition. Note 
that the c and d coefficients are consistently used in connection 
with the dT/dp and dV/dq term respectively. The computa¬ 
tional cost of Sis is approximately proportional to the number 
of derivative evaluations (stages). The number of stages is de¬ 
fined as s = k — 1 for symmetric coefficients due to the first 
same as last property, which implies that the first derivative 
evaluation is the same as the last derivative evaluation of the 
previous iteration. 


by the entries Q;, ^ in Table 2 where h denotes the intermedi¬ 
ate step number. This decomposition follows directly from Eq. 
6. A similar decomposition can be constructed for the ABA 
composition. The position contribution £a for each order A is 
introduced as: 

Cat a = Q s , A T A ,A = 0,l,...,A max (7) 

where s is the number of stages. The new position can then be 
calculated as the sum of all position contributions: 


Table 1: ABA and BAB composition for symmetric coefficients 
ABA BAB 


di = <4_i 

Cl c k 

di = d k 

Cl = c fc ! 

d 2 = <4- 2 

c 2 = Cfc—1 

di = <4-1 

c 2 = Ct_ 2 

d k _ 2 



Cfc—2 

dk- 1 

Cjt-1 

dk- 1 

Cjfc-1 

d k = 0 

Cfc 

dk 

Ck =0 


III. Methods 

This section first describes a decomposition method which 
follows directly from the composition of the simple harmonic 
oscillator. It is then shown that this decomposition can be 
used to satisfy general fourth order conditions and higher 
order conditions for the simple harmonic oscillator. 


7 1 ma x , >. 

% = E t a ^ a + o (t a ”“ +1 J 


Suppose that a known function q(t) exists such that q(to) = q a 
and q(to + r) = q b + o( r Am «* +1 ) then q b is expressed by the 
power series: 


<7 ( f o + T ) 


lW T l I c l"( t o)^2 I , q (Amaxl ( f o)^A, 


= <7« + + -4t^T z + ... + 




= Co + Ch 1 + £ 2 t 2 + ... + C\ max T Amax + 0 


0 


If the SI is of Nth order then 


<? (A) (to) 


t a = £ a t a ,A = 0,1 


for any function q(t) that qualifies as a time evolution in 
H(p,q) = T(p) + V(q). The simple harmonic oscillator (SHO): 

H = (11) 


A single iteration from (q a , p a ) to (q b , p b ) using the s-stage 
BAB-composition (Eq. 4) can be expanded into: 

pi =p a -Td 1 ^(q a ) 
q 1 = q a + TC^(p 1 ) 

V 2 = pi-Td 2 ^r (qi) 

(5) 

<?s = (? 5 _1 + TC S |^ (p s ) = q b 
p s+1 = p s - Td s+ 1 ^ (q s ) = p b 

Suppose now that dT ( p)/dp = p/m which applies to many 
classical systems including the simple harmonic oscillator. Eq. 
5 then yields: 

q 1 = q a + rci^ + (q a ) 

<72 = < 7 i (1 + C 2 ) - C 2 q a + t 2 D 2 ;;A 1 |^ (qffi 
<73 = <72 (1 + C 3 ) - C 3171 + r 2 D 3 m _1 |^ (q 2 ) 

( 6 ) 

<?s = < 7 s 1 (1 + C S ) - C s q s -2 + T-Dsitr 1 ^ ( <7s _ 1 ) 

Ps +1 = - Td s+i% (< 7 s) 

where C„ = c n+ i/c n and D„ = —c n d n ■ The first same as last 
property allows substitution of p a by p s+ i yielding an expres¬ 
sion that is independent of p. Eq. 6 yields significantly larger 
truncation errors than Eq. 5 and is therefore not suitable for 
practical implementation. However, Eq. 6 allows separation of 
t a terms for any order A = 0,1,..., A max where A max is the max¬ 
imum number of error terms which appear by decomposition. 
For instance consider the three stage decomposition expressed 


is enforced by applying the condition Qi h \ = dV/dq(Qu,\) in 
Table 2. The SHO yields the well-known analytical time evolu¬ 
tion of q: 



Using Eq. 10 where £o,£i,..., £6 are evaluated in Table 2 the 
following identities must be true for symmetric fourth order 
splitting coefficients: 

h = Co = <7« 

v.\fl* = h = p 4 ( 2 ci+c 2 ) 

_ !! fl = £2 = “IF (til +d 2 ) ( 2 ci +C 2 ) (13) 

- A! (!) = £ 3 = -2^2^ (Cl + c 2 ) 

jr (!) a = £4 = Cid 2 ^//f ( 2 cidi + lc 2 di + c 2 d 2 ) 

Elimination of the constants a and b yields the identities: 

1 = 2 (di + d-2)(2ci + C 2 ) 

1 — \2c\d-2(c\ +C2)/(2c;l + C 2 ) (14) 

1 — 24 cid 2 ( 2 cidi 2 c 2 d\ H- ^2^2) 

It is readily found that the only real solution satisfying Eq. 14 
is the one originally found by [1]: 

di = d 4 = 2 f 2 ! 2 l 73 \ =\-d 2 = \-d 2 

1 ’ 1 /, x ( 15 ) 

Cl = c 3 = 2^5173 = 2 (! - c 2 


A closer examination of Table 2 reveals that m, k and q a only 
appear as scalars of £q,£i, £a„,„ the r°, and that t 2 and r 4 
separations can be expressed as in Table 3. 
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Table 2: BAB decomposition for symmetric coefficients where Q, = C/, + 1, Qh A = m 1 dV/dq(Qh,\) and Qh,x denotes the table value at indices (h, A) 



T° 

T 1 

T 2 

T 3 

T 4 

T 5 

T 6 

0 

1 

2 

3 

Qo,o 

Ql,0 C 2 - C 2 Qo,0 
Ql,oC 3 “ C 3 Qi,o 

Cipam^ 1 

Ql,\C 2 

Q24C3 — C3Q14 

DiQo.o 

Ql,2C2 + D2Q1,0 

Q 22 C 3 — C 3 Q 12 + D 3 Q 2/ 0 

D 2 Q 1 A 

Q23C3 + D3Q2T 

DlQl,2 

Q24C3 + D 3 Q 2/ 2 

D 3Q2,3 

D 3 Q 2 A 


Table 3: Decomposition of even orders expressed in terms of the splitting 
coefficients, m and q a are set to 1 


h A 

T° 

T 2 

T 4 

0 

1 



1 

1 

1 1 

- E Ci E dj 

i= 1 /=! 


2 

1 

2 i 

- E c; E dj 
i= 1 y=i 

1 1 
— Y2 Qm,2^m+1 YL C k 
m= 1 k=m+ 1 

h 

1 

h i 

- E Ci E dj 
«'= 1 M 

h -1 h -1 

— Yh Qm,2^m+1 YY Ck 

m =1 k=m +1 


Symmetric coefficients omit odd order error terms [2], im¬ 
plying that the following identities must be true to satisfy the 
fourth conditions for the SHO: 


£2 

u 


<7"(0) _ 1 

2 ! “ 2 ! 

<f""(0) _ 1 
4! “ 4! ' 


h i 

- E Ci E dj 
»■=1 y=i 

ft-1 

E Qm,2“m+1 

m=l 




E c* 

/c=m+l 


( 16 ) 


It has previously been shown that Eq. 16 is in fact general 
second and fourth order criteria [6]. This implies that split¬ 
ting coefficients satisfying Eq. 10 to at least fourth order also 
satisfy Eq. 2 to at least fourth order. Compositions with more 
stages and higher order SHO criteria can be exploited for sig¬ 
nificantly reducing higher order error terms. In this sense, 
the higher order SHO criteria are used only as a constraint 
for obtaining efficient splitting coefficients for more general 
Hamiltonian systems. 


Analytical evaluation of splitting coefficients for more than 
three stages quickly becomes hopeless. Instead a numerical 
routine has been developed to satisfy the multi-objective op¬ 
timization problem with variables c and d and minimization 
objective k: 

k a = \qW (U) - \ = 0,A = 1,2,...,A H (17) 

where Ah is referred to as the SHO order constraint and k a 
for A = 1,2,..., A h are the minimization objectives. Numerical 
experiments suggest that for A h > 5 every additional stage 
yields two increments of A max and a single increment of A h 
for which at least one real solution exists. 


The BAB decomposition is numerically evaluated with the 
choice of constants in Eq. 12: a = — b = 1 and m = k = 1 
corresponding to the initial conditions q a = — p a = 1. The 


objectives k 3 and xy are always 0 if E/Li c i = E?=i dj = 1. 
Additionally these identities eliminate one c and d variable 
from the minimization procedure. The minimization is per¬ 
formed using an adaptive simplex method [7] minimizing 
only one objective K max = max{K\,K 2 ,...,K A ). Solutions are 
provided in the following section with minimization criteria 
K max < le — 128 calculated using software-implemented ar¬ 
bitrary precision arithmetic (2048bit, 256bit exponent). Four 
unique sets of splitting coefficients are presented in the fol¬ 
lowing section. The non-unique solutions are picked among 
thousands of solutions based on the sum of the absolute val¬ 
ues of the splitting coefficients and higher order SHO errors. 


IV. Results 

A. Obtained splitting coefficients 

Splitting coefficients with five to nine stages are presented 
in Table 4 in array format for convenient implementation for 
various programming languages. 

Two different kinds of BAB solutions are evaluated: 

• BAB: Solution to Eq. 10 using BAB decomposition (see 
Table 2). 

• BAB': Solution to Eq. 10 using ABA decomposition with 
coefficients cq = c s+ i = 0 in Eq. 3. 

The difference between the BAB and BAB' solutions is re¬ 
vealed by Eq. 5 where the error of p s+ \ remains unaddressed 
for the BAB solution as opposed to BAB' solution where the 
error of q s+ \ remains unaddressed. Thus the BAB and BAB' 
solutions possess different constraints which may influence 
their performance in either direction. 

BAB and BAB' solutions should be implemented accord¬ 
ing to Eq. 4. Unique solutions satisfy both the ABA and BAB 
decompositions and can therefore be implemented as either, 
although the performance of the presented solutions depends 
on the applied scheme. 

The naming convention of the obtained methods is 
[scheme]s[stages]o[order]H, for instance BAB's9o7H for a BAB- 
optimized scheme with nine stages and seventh harmonic or¬ 
der. 
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B. Numerical tests 


Common methods used to benchmark Sis for integration 
of autonomous Hamiltonian systems include [4] [5] [8]: 

• The maximum Hamiltonian error: max(\AH(t) /Hg|) 

• The mean Hamiltonian error: mean(\AH(t)/Hg\) 

• The accumulated error of integrals 


where Hg is the initial Hamiltonian. Each of these benchmark 
methods are applied throughout this section. 


The following 10 symmetric Sis are benchmarked for three 
autonomous Hamiltonian systems: The simple harmonic os¬ 
cillator, the Henon Heiles system and the Kepler problem: 


• Ruth [ABA] 3-stage, 4th order [1] 

• s5odr4 [ABA] 5-stage, 4th order [9] 

• ABA104 [ABA] 7-stage, 10-4 generalized order [5] 

• ABA864 [ABA] 7-stage, 8-6-4 generalized order [5] 

• ABA1064 [ABA] 8-stage, 10-6-4 generalized order [5] 

• Yosh s7o6 A [ABA] 7-stage, 6th order [2] 

• ABAs5o6H [ABA] 5-stage, 4th order 

• BABs7o7H [BAB] 7-stage, 4th order 

• BAB's8o7H [BAB] 7-stage, 4th order 

• BAB's9o7H [BAB] 9-stage, 4th order 


where [ABA] and [BAB] denote suggested ABA- and BAB- 
implementation respectively. It should be noted that the 
generalized order schemes are designed for perturbed Hamil¬ 
tonian systems H = Ha + eHg with e <C 1. It should also 
be noted Yosh s7o6 A yields the smallest errors of the three 
stage-optimized sixth order integrators presented by [2], and 
that this integrator is here evaluated to reference higher order 
accuracy and stability. 


Numerical tests are performed using double precision 
floating point format with compensated summation. Com¬ 
pensated summation drastically reduces truncation errors 
for long-term symplectic integration (see [9] and references 
therein). All numerical integrations are performed using Zym- 
plectic (v.1.01.00) from Zymplectic.com, which is a numerical 
framework designed for numerical integration and benchmark 
of separable Hamiltonian systems. 
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Figure 1: Integrator benchmark using the SHO in the time interval t 6 
[0,500]. Top: Maximum Hamiltonian error. Bottom: Average 
relative Hamiltonian error. The dashed lines show benchmark pro¬ 
files for integrators of corresponding color using extended preci¬ 
sion. The horizontal axes correspond to r/s. All axes are in base 
10 logarithmic units. Legend applies to both graphs 


The simple harmonic oscillator 

Consider the Hamiltonian of the one-dimensional SHO (Eq. 
11) with k = m = 1 and initial conditions (p,q) = (0,1). 
Figure 1 shows a benchmark profile of the 10 integrators for 
the SHO. Every point on each line corresponds to a complete 
integration in a given time interval. The slopes of the lines 
read the integrator order. It is noticable that the integrators 
ABAs5o6H A, BABs7o7H, BAB's8o7H and BAB's9o7H be¬ 
have as sixth order integrators. 


The ABA- and BAB-decomposition from the previous sec¬ 
tion can be utilized to determine error contributions of differ¬ 
ent orders. More precisely, K\ can be evaluated for all A to 
resolve the magnitude of higher order errors for the SHO. This 
can be used to characterize higher order error terms which of¬ 
ten diverge rapidly for high order stage-optimized integrators. 
While these error terms strictly apply to the SHO, they still 
cover general equations present in the expansion of Eq. 2. 
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Figure 2: SHO decomposition of selected integrators for both the ABA- (cir¬ 
cle marker) and BAB-composition (solid line). Lines and markers 
of the same color intersect exactly at even orders 


Figure 2 shows position error terms of selected integrators for 
the SHO. Notice that the obtained BAB methods have more 
error terms when implemented as ABA methods. Further¬ 
more the seventh harmonic order of BAB's8o7H eliminates 
the seventh order error term of the position coordinate for the 
BAB implementation as opposed to the sixth harmonic order 
integrator ABAs5o6H A. The largest integer value of A with 
no integration error reads the SHO integration order. Stage- 
optimized integrators of higher orders generally have rapidly 
diverging higher order error contributions. This implies that 
these integrators are not suitable for integrating the SHO with 
large time steps. 


The Henon Heiles system 

The Henon-Heiles system is one of the most studied Hamilto¬ 
nian systems due to its applications in celestial mechanics, its 
non-integrable properties, chaotic dynamics and fractal struc¬ 
tures. Consider the Hamiltonian of the Henon-Heiles system 
with initial conditions [q x , q y , p x , Py\ = [0.3,0.0,0.0,0.4]: 


H = 


vl+v] qf+qf 

2 H 5— 


+ <7 x 2 q y - \q y 3 


(18) 


where in this case H = 1/8, inside the chaotic domain. Fig¬ 
ure 3 shows a benchmark profile of the 10 integrators for the 
Henon Heiles system. The Henon Heiles system can be consid¬ 
ered as a perturbed SHO for small values of q x and qy. For this 
reason it appears that especially BAB's8o7H and BAB's9o7H 
exhibit a sixth order descent of error for large r. 



time per evaluation 



time per evaluation 


Figure 3: Integrator benchmark using the Henon Heiles system in the time 
interval t € [0,500]. Top: Maximum Hamiltonian error. Bot¬ 
tom: Average relative Hamiltonian error. The horizontal axes cor¬ 
respond to t/s. All axes are in base 10 logarithmic units. Legend 
applies to both graphs 


The Kepler problem 

Symplectic integration has been widely used for high accu¬ 
racy simulations of the Solar System (see [8] and references 
therein). It has been proposed that Mercury due to its fast or¬ 
bital period and high eccentricity is limiting the performance 
of Sis when integrating the Solar System [8]. Consider the 
Hamiltonian of the planar Sun-Mercury system: 

zj __ IIpII 2 nm s m M ~ l|p|| 2 nm s m M ncn 

H ^iqT ~ ^ _ G_ Hqir [ ’ 
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where M is the reduced mass, G is the gravitational constant, 
and mg and m m are the masses of the Sun and Mercury re¬ 
spectively. The planar Sun-Mercury system is here evaluated 
using data from nssdc.gsfc.nasa.gov where the aphelion is cho¬ 
sen as the initial Sun-Mercury distance yielding the minimum 
orbital velocity. 

Figure 5 shows a benchmark profile of the 10 integrators 
for the Sun-Mercury system. It is noticable that ABA104 
shows the best performance in regard to the maximum error 
and BAB's9o7H shows the best performance in regard to the 
average error. 

The accumulated error in a two body system (Sun-Mercury) 
manifests as a secular advance of the orbital ellipse similar to 
the precession of perihelion induced by general relativity or 
the gravitational pull of other planets. The mathematical na¬ 
ture of the secular SI error is described by [10]. The secular 
advance of the orbital ellipse can be evaluated through the 
rotation of the Laplace-Runge-Lenz vector A: 



time per evaluation 


A = p x L — r 


( 20 ) 


where r = r/r, L = rxp, pis the momentum vector and r is 
the position vector with length r. 
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Figure 4: Advance of perihelion (absolute rotation velocity) as a function of 
the number of evaluations per orbit. Linear fit 95% confidence in¬ 
tervals of the perihelion advance are significantly smaller than the 
dot size. The integrators ABAs5o6H, ABA104 and BAB's8o7H 
achieve approximately the same performance, obfuscating their 
plots 



Figure 5: Integrator benchmark using the planar Sun-Mercury system in 
the time interval t 6 [0,10yr]. Top: Maximum Hamiltonian er¬ 
ror. Bottom: Average relative Hamiltonian error. The horizontal 
axes correspond to t/s in units of hours. All axes are in base 10 
logarithmic units. Legend applies to both graphs 


The orbital precession advances linearly with time t [10]. The 
rotation 8 is introduced as the angle of A implying that 0(f) 
is a linear function and is chosen such that 0(0) = 0. Fig¬ 
ure 5 shows the perihelion precession d8/dt which has been 
obtained numerically through linear regression of 0(f) over 
numerous orbital periods. Comparison of Figure 4 and 5 sug¬ 
gests that the accumulated error is in better agreement with 
the average error than the maximum error. 
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Figure 6: Phase space illustration of the substeps performed over two iterations by Yosh s7o6 A (A)(C) and BABs7o7H (B)(D). The integrators are applied to 
the SHO (A)(B) with initial conditions ( q,p ) = (0,1) and time step t = n/i, and the Henon-Heiles system (C)(D) in the y-plane ( q x = p x = 0) 
with initial conditions ( q,p ) = (qy,Py) = (0.4,0.4) and time step t = 2. Blue circles and Xs indicate the initial and final positions respectively 
for each iteration. Blue dots indicate substeps with positive c; and df coefficient for calculating the next substep. Red dots and green dots indicate 
substeps with negative c; and dj coefficients respectively for calculating the next substep. A green dot within a red dot indicates that both c ; and 
di are negative. The contour colors display the Hamiltonian with values given by the colorbars. (C) shows an unstable integration where H is not 
conserved. (A) (B) and (D) show stable integrations which preserve H indefinitely 


V. Discussion 

A thorough benchmark analysis of 10 different integrators 
has been performed for the SHO, the Henon Heiles system 
and the Kepler problem. All numerical tests suggest that the 
obtained integrators in spite of an increased number of stages 
greatly and consistently improve the integration accuracy and 
stability. In fact, the new schemes outperforms the acclaimed 
stage-optimized integrator Ruth by several orders of magni¬ 
tude. The higher order behavior for SHO-perturbed systems 
is an additional perk of the obtained integrators. This prop¬ 
erty may be particularly useful in fields of molecular dynam¬ 
ics where symplectic correctors may not be suitable. Opti¬ 
mized fourth order schemes appear to be useful for achiev¬ 
ing high numerical stability allowing large time steps with 
minimal integration error. High stability is particularly im¬ 
portant for systems where the complexity of H varies with 
time: Close encounters in the gravitational N-body system in¬ 


evitably cause high errors ultimately waiving symplectic prop¬ 
erties of Sis. Higher order integrators are generally more con¬ 
strained and usually involve several negative splitting coeffi¬ 
cients with large absolute values. The performance breakpoint 
of higher order integrators can be visualized by considering 
the executed substeps. Figure 6 illustrates the substeps per¬ 
formed by two seven stage Sis (Yosh s7o6 A and BABs7o7H) 
for the SHO and Henon-Heiles system. Figure 6 reveals that 
Yosh s7o6 A (A) (C) performs large forward and backward 
steps compared to BABs7o7H (B) (D). The steps are too large 
in (C) for the Hamiltonian to be preserved. BABs7o7H and 
other splitting coefficients obtained in this paper generally per¬ 
form small near-forward steps. Note that some of the methods 
presented in Table 4 only have small coefficients for either d or 
c. For this reason the integrator performance may change dras¬ 
tically if the roles of T and V are switched. The performance of 
the obtained integrators generally declines for systems where 
T(p) deviates significantly from a squared kinetic energy. 
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VI. Conclusion 


Efficient fourth order splitting coefficients have been ob¬ 
tained by satisfying higher order conditions for the simple 
harmonic oscillator. Numerical tests indicate that the obtained 
integrators outperform the three stage fourth order integrator 
by R. Ruth as well as the more optimized integrators by W. 
Kahan and S. Blanes. The obtained integrators also compare 
favorably with higher order integrators by H. Yoshida for 
moderate time steps. The integrators are particularly suitable 
for particle systems that require high numerical stability and 
accuracy at large time steps. The best performance is gen¬ 
erally achieved by the methods BAB's9o7H, BAB's8o7H and 
ABAs5o6H A. 


It has been found that harmonic order conditions higher 
than seven often introduce more negative splitting coefficients. 
A larger number of stages presents more degrees of freedom 
resulting in more inappropriate solutions. Additional con¬ 
straints should be pursued to obtain efficient splitting coef¬ 
ficients with more than nine stages. 

The presented decomposition can be used to satisfy gen¬ 
eral fourth order criteria and higher order conditions specif¬ 
ically for the simple harmonic oscillator. This work suggests 
that general purpose high accuracy, high stability explicit sym- 
plectic integrators can be pursued by approaching multiple 
higher order conditions. This has been achieved here by satis¬ 
fying individual high order entries. 


References 

[1] Fourth-order symplectic integration. Forest, Etienne and 
Ruth, Ronald D. 1, 1990, Physics D, Vol. 43, pp. 105-117. 

[2] Construction of higher order symplectic integrators. 
Yoshida, Haruo. 5-7, 1990, Physics Letters A, Vol. 150, 

pp. 262-268. 

[3] Forward Symplectic Integrators for Solving Gravita¬ 
tional Few-Body Problems. Chin, Siu A. and Chen, C. 
R. 3-4, 2005, Celestial Mechanics and Dynamical Astron¬ 
omy, Vol. 91, pp. 301-322. 

[4] Wisdom, Holman, M. and Touma, J. Symplectic 
Correctors, [book auth.] Jerrold E. Marsden, George 
W. Patrick and William F. Shadwick. Integration Algo¬ 
rithms and Classical Mechanics, s.l. : Providence, RI, 
1996, Vol. 10, pp. 217-244. 

[5] New families of symplectic splitting methods for numer¬ 
ical integration in dynamical astronomy. Blanes, Sergio, 
et al. June 2013, Applied Numerical Mathematics, Vol. 
68, pp. 58-72. 


[6] Optimized fifth order symplectic integrators for orbital 
problems. Tselios, Kostas and Simos, T.E. 1, 2013, Re¬ 
vista Mexicana de Astronomia y Astrofisica, Vol. 49, pp. 
11-24. 

[7] Implementing the Nelder-Mead simplex algorithm with 
adaptive parameters. Gao, Fuchang and Han, Lixing. 
1, 2012, Computational Optimization and Applications, 
Vol. 51, pp. 259-277. 

[8] High precision symplectic integrators for the Solar Sys¬ 
tem. Farres, Ariadna, et al. 2, 2013, Celestial Mechanics 
and Dynamical Astronomy, Vol. 116, pp. 141-174. 

[9] Composition constants for raising the orders of uncon¬ 
ventional schemes for ordinary differential equations. 
Kahan, William and Li, Ren-Cang. 219, 1997, Mathemat¬ 
ics of Computation, Vol. 66, pp. 1089-1099. 

[10] Physics of symplectic integrators: Perihelion advances 
and symplectic corrector algorithms. Chin, Siu A. 3, 
2007, Physics review E, Vol. 75. 


8 



Table 4: Obtained splitting coefficients with 77 digits. The methods are organized according to their stage count. 


Unique ABAs5o6H, A, B and C 


d(1) = 0.1558593591762168313166117535752091422239663993391011462498104831549442591694 
d(2) =-0.0070254990919573173514483364758218294773716640092220571342056284758867609611 
c(1) =-0.6859195549562166768601873150414759494319985863677163820719179393682014399373 
c(2) = 0.9966295909529363159571451429325843698583459772292551181721475637244006507927 

d(l) = 0.4020196038964999834667409950496227775945673320979099323902806525851620445492 
d(2) = 0.5329396856308538150258772262086702929451721575835842834460326556965220312130 
c(1) = 0.9110842375676615218574607388486783304139753525628699898390474132061253024968 
c(2) = 0.1740059542332660799009374186088931171982348451547482386207462271424421679090 

d(l) = 0.1868565631155112597511173758337610451623768791420295598869906080256347098408 
d(2) = 0.5520581660514781484261043096825685955052553493857487316732455515112095793516 
c(1) = 0.5642486163110637621453746447826190031465518453448216443978248524452914255263 
c(2) =-0.2393627021773294286793711975145735718917010075899623225091609656425715483488 


d(3)=0.5-d(l)-d(2) ; d (4) =d (3) ; d (5) =d (2) ; d (6) =d (1) ; 
c (3) =1-2* (c (1) +c (2) ); c (4) =c (2); c (5) =c (1); 

Unique BABs6o7H (more exist), BABs6o5H and BAB's6o5H 


d(l) = 0.0832701092493097690276300822599156817795619881080575430174826369500044839553 
d (2) = 0.3997273690963360211284395 9200077 955505750605316347 937480202072 8897 63444 394 68 
d(3) =-0.0541842778124726964199287659702152862181671805554302053695494244408226729818 
c(1) = 0.2475471587650765967910125296669232190787926795528258860075742877866898482465 
c(2) = 0.5446579217808193419580029125986805136192611468678745304306198457355253495088 

d (1) = 0.06588315331611550217 943712 97 62 994 92142112706414503678821086524542252382 92357 
d(2) =-0.6711629060948253965117521242801468651670183829736696004743743104248379033034 
d (3) = 0.97367031007253504 98414312 65155085719113121893230878832080676206300409632 0895 
c(1) = 0.2265023974336291596186923088995152371194987043433278784212774210853229477086 
c(2) =-0.0047799986678794678665602622568725658855054645768977416258774175978957628102 

d(l) = 0.0650508268637574949487516678539036744380576078003520231297474895927182047842 
d(2) =-0.3948051939117155639582651907195511796839512131373933326629623631591978696178 
d(3) = 0.6918498547904058960782554213200966000604457266088855079657967408955821538731 
c(l) = 0.2328962665845291347812910553597276545034489573682700034501734659308763580659 
c(2) =-0.0111617638003721094728940473306267483522869816097800738075975236192041340802 


d (4) =1-2* (d (1) +d (2) +d (3) ) ; d (5) =d (3) ; d (6) =d (2) ; d (7) =d (1) ; 
c(3)=0.5-c(l)-c(2); c (4) =c (3); c (5) =c (2) ;c(6)=c(l); 

BABs7o7H and BAB's7o6H 


d(1) = 0.0638745574250616045658401356462756092272737349204789877616691621039130680037 
d(2) =-0.0650239777505938311516598494765811300128929849501107531440553736739769298894 
d(3) = 0.2509446105745547370613575645855473357282136355718617210088090794709222342775 
c (1) = 0.275278172 905 9777393394 97871044869078212521501894 918608507532560534 852 6197 7 56 
c(2) =-0.0843138705589167473554015820986490036832890668438279781819362930106920807542 
c(3) = 0.1674497222006475614401177016323447087805836086414469568091358611098423440220 

d(1) = 0.0522155297747848201407012160969040693245471580104797248381281194964273517726 
d (2) =-0.0824 97255852 95 6141213191193771742051416272833968105 6503508469680313691406287 
d(3) = 0.3285541797987193353601113204079269672646845923663727576986276602617960257026 
c(l) = 0.2487563308365098625528031803769571289196558939258433219240690943143969198069 
c(2) =-0.0651011247076581799932061212576878177123945470202647856511921757791017775052 
c(3) = 0.2480624780675545152650672751613106579864581926645260137078906816928505888862 


d (4) =0.5-d (1) -d (2) -d (3) ; d (5) =d (4) ; d (6) =d (3) ; d (7) =d (2) ; d (8) =d (1) ; 
c (4) =1-2* (c (1) +c (2) +c (3) ) ; c (5) =c (3) ; c (6) =c (2) ; c (7) =c (1) ; 

BAB's8o7H 


d(1) = 0.0538184115480034769403763798524605188562842390760879592632218376015166638395 
d (2) = 0.164874332 691047236101480908531705 94252 99121141031052090901977 95251398487 8990 
d(3) = 0.3895399407808198068744134256203146340834631254960864069726823667050364522355 
d(4) =-0.2288957415563594299572505173565338312542463595622272333825061768110435645417 
c(1) = 0.1486140577445185629163082471176700173109512976367237631150576219945233462284 
c(2) = 0.1071986675806227950500566279939336794589433458464489776124879870581484936262 
c(3) =-0.0149646736494517061945681450558142918831874360003431672632178480031079700216 


d (5) =1-2* (d (1) +d (2) +d (3) +d (4) ) ; d (6) =d (4) ; d (7) =d (3) ; d (8) =d (2) ; d (9) =d (1) ; 
c(4)=0.5-c(l)-c(2)-c(3) ;c (5) =c (4) ;c (6) =c (3) ;c(7)=c(2) ;c(8)=c(l) ; 


BAB's9o7H 


d(1) = 0.0464929004396589154281717058427105561306160230440930588914036807441235817244 
d(2) = 0.1549010127028879927850680477816652638346460615901974901213193690401204696252 
d (3) = 0.3197054828735 917137 611074311771339117 602 99488424509122033334 0037 84161608504 8 
d(4) =-0.1929200088157132136865513532391282410293753210475133631464188500663304857888 
c(1) = 0.1289555065927298176557065467802633438775379080212831185779306825670371511433 
c(2) = 0.1090764298548827040268039227200943338187149719339317536310302288046641781422 
c(3) =-0.0138860356804715144111581981849964201100030653749527555344377031679795959892 
c(4) = 0.1837549745641803566768357217228586277331494085368674804908537743649129597425 


d (5) =0.5-d (1) -d (2) -d (3) -d (4) ;d (6) =d (5) ;d (7) =d (4) ;d (8) =d (3) ;d (9) =d (2) ;d(10) =d (1) ; 
c (5) =1-2* (c (1) +c (2) +c (3) +c (4) ) ; c (6) =c (4) ; c (7) =c (3) ;c (8) =c (2) ;c (9) =c (1) ; 
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